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Abstract 

A discontinuous Galerkin method has been developed for strain gradient- 
dependent damage. The strength of this method lies in the fact that it allows 
the use of C° interpolation functions for continuum theories involving higher- 
order derivatives, while in a conventional framework at least C interpolations 
are required. The discontinuous Galerkin formulation thereby offers significant 
potential for engineering computations with strain gradient-dependent models. 
When using basis functions with a low degree of continuity, jump conditions 
arise at element edges which are incorporated in the weak form. In addition 
to the formulation itself, a detailed study of the convergence properties of the 
method for various element types is presented, an error analysis is undertaken, 
and the method is also shown to work in two dimensions. 



1 Introduction 



The nucleation and growth of cracks can be described by continuum damage 
mechanics. In this approach, an internal variable is included in the cons t itutiv e 
law to represent the evolution of microstructural damage l|Kachanovt Il958j) . 
The material moduli degrade with the increase of this parameter, which could 
be either scalar or tensorial in character. A scalar internal variable is commonly 
used to model damage degradation for cases in which the material remains 
isotropic. A tensorial form of the internal variable is used for the anisotropic 
case. 

Damage degradation can manifest itself in progressive material softening, for 
which reason numerical results based upon classi cal cont i nuum mechanics ar e 
characterised by a pathological mesh dependence llWillaml . Il98l iBazantl Il98r3) . 
As the material softens, deformations localise within bands o f finite thickness , 
which is viewed as the smearing out of a crack (|Rashidl Il968l iBazant and Old . 
119831) . However, the width of bands is intimately related to the element size. 
As the mesh is refined the band width decreases, tending to zero in the limit. 
The result of this is seen in load-displacement response with slopes that are 
increasingly negative with mesh refinement in the softening regime. The origin 
of this pathology lies in the loss of ellipticity of the material tangent modulus 
tensor for softening inelastic continuum models in the absence of rate effects 

The last two decades have witnessed a great deal of work in 'regularised' con- 
tinuum models to avoid the loss of ellipticity in the presence of material soften- 
ing. Among these are the strain gradient models llAifantislll984t IColeman and Hodgdonl 
ll985tlTriantafvllidis and Aifantislll986t|P eerlings et allll996|). nonlocal model s 
l|Bazant and Piiaudier-Cabotl ll98St Ide Vree et all Il995t iBorino et all l2003lh 
locali sation limiters l)LasrT^mdFieTvtcnko|.ll998j) . and Cosserat models l|de Borst and Sluvsl 
199l|, to identify but a fraction of a vast body of literature. 

By various techniques, each of these models introduces an intrinsic length 
scale to the continuum theory. When the material softens, the localisation band 
width is controlled by this length scale, rather than the element size. As a 
result, the mesh size-dependency, described above, is eliminated. 

In this paper we aim to address numerical issues associated with some strain 
gradient models. This class of models involves strains and strain gradients as 
kinematic terms. In some models work conjugate couple stresses are also iden- 
tified. Dimensional considerations lead to the introduction of intrinsic length 
scales associated with the strain gradients. The presence of strain gradient 
terms alters the mathematical character of the governing differential equations, 
elevating them to a higher order (at least locally). A numerical complication 
arises from the higher order character of the governing differential equations, as 
standard C° interpolation functions often do not furnish sufficient continuity. 

The obvious approach, that is to employ interpolation functions with higher- 
order continuity (C 1 and higher), has not proven fruitful. Finite element meth- 
ods with these interpolations are expensive and difficult to construct. There 
exists some work with mixed formulations for this purpose l|Shu et all Il999t 
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Zervos et all 12001). However, the very simplest two-dimensional elements in- 
volve up to 28 degrees of freedom, placing them beyond the realm of prac- 
ticability. Some authors have applied element-free Galerkin met hods, on the 
basis of the arbitrary degree of continuity that they can provide ijAskes et all 
2000). However, this class of methods introduces complications in the applica- 
tion of boundary conditions and is somewhat less efficient than the finite element 
method, and its implementation involves a reformulation of the dominant finite 
element architecture, making its widespread use less attractive. In addition, 
just as high-order continuity can be difficult to achieve, excessive continuity can 
also be a drawback for many models, particularly at internal boundaries where 
the order of the governing differential equation changes. 

In this work we seek to expand upon recent developm ents in discontinu- 
ous Galerkin meth ods for strain gradient-dependent theories llEngel et al 
IWells et aD . 120041) . The great advantage of this class of methods is the ability 
to use C° interpolation functions for the displacement when solving continuum 
theories involving higher-order derivatives. In cases where other considerations 
make it advantageous to also interpolate strain-like quantities, the approach can 
allow interpolated fields that are discontinuous across elements ll Wells et all 
120041) . 

The content of this paper is organised as follows. The strain gradient damage 
model of interest is described in Section O with special attention paid to the 
boundary conditions. The Galerkin formulation for this model is presented in 
Section|3| The convergence behaviour of the model in one dimension is examined 
both analytically and numerically in Section^] and the section is concluded with 
one- and two-dimensional examples. Conclusions are then drawn in Sectional 



2 Gradient damage 

Consider a body fl, which is an open subset in M. n , where n is the spatial 
dimension. The boundary of the body is denoted by T = <9fi, and the outward 
unit normal to T is denoted by n. We consider quasi-static equilibrium of the 
body, which is governed by 

V • <r + / = in 0, (1) 

an = h on Th, (2) 

u = g onr g , (3) 

where er is the stress tensor, / is the body force, h is the prescribed traction 
on the boundary Th and g is the prescribed displacement on the boundary Y g . 
The boundary T is partitioned such that T — T g U Th and T g n I\ = 0. For 
isotropic damage, it is sufficient to introduce a single scalar damage variable, 
D, satisfying < D < 1. Considering that D = implies no damage and D = 1 
corresponds to a completely damaged material point, the stress-strain relation 
is of the form 

o-=(l-D)C:e, (4) 
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where C is the elasticity tensor and e is the strain tensor, e = V s it. The 
scalar damage variable, D, is related to a scalar strain measure, e, through the 
Kuhn- Tucker conditions: 

e - k < 0, k > 0, k(e - k) = 0, (5) 

where -D = D(n). The scalar strain measure e of interest here arises from a 
so-called explicit gradient damage formulation. In this approach, the strain 
measure is defined to be 

e = £cq + c 2 V 2 e oq , (6) 

where c is an intrinsic length scale, and e cq is a suitably chosen invariant of the 
local strain tensor. 

For material points at which damage is developing, the governing equation 
is nonlinear and fourth-order in the displacement. In undamaged or unloading 
regions (k = 0), the governing equation is linear and second-order. Recognising 
the fourth-order form of the differential equation (at least in sub-regions of f2) 
leads to the question of appropriate boundary conditions. Possible boundary 
conditions (in addition to the usual boundary conditions on the displacement 
and traction) include: 

£ cq = q on r £ , (7) 

\7e cq -n d = r onr E -, (8) 

(V (Ve oq ) n d ) -n d = s on T e „ , (9) 

where r e U IV U T e » = T d , T E n IV = 0, I\ H T e >, = 0, T E , n T e » = 0, and n d 
indicates the unit outward normal to the boundary of the damaged domain, IV 
These boundary conditions supplement the standard displacement and traction 
boundary conditions on T. If the entire body is damaging, then T d — T. 
Commonly however, damage is localised. In this case interface conditions arise 
on the boundaries which are internal to O (the boundaries to sub-domains where 
k > and k — meet). When the body force is smooth, continuity of u and 
crn d are natural continuity conditions. At the interface between damaging and 
elastic regions, an additional condition is required on the damage domain. The 
need for extra boundary conditions is made particularly clear by the analytical 
solution presented in later in Section 14.2.11 In the adopted formulation, this 
extra condition will come from implied continuity of X?e eq -n d across the damage- 
elastic boundary. 

3 Galerkin formulation 

In developing a weak formulation for eventual finite clement solution, the equi- 
librium equation Q and the equation for e © are considered separately. The 
nonlinear fourth-order equation resulting from insertion of the constitutive equa- 
tions into the equilibrium equation could potentially be cast in a weak form, 
with integration by parts applied twice. The formulation would be specific 
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to the chosen dependence of damage on k, which is typically complex. Fur- 
thermore, the complexities associated with a moving elastic-damage boundary 
(which is the interface between second- and fourth-order sub-domains) would 
be significant. Hence, it is convenient to treat the two equations separately. 
The body £1 is partitioned into n c \ non-overlapping elements f2 e such that 



n = (J n e , (io) 



e=l 



where fi e is a closed set (i.e., it includes the boundary of the element). The 
elements Q e (which are open sets) satisfy the standard requirements for a finite 
element partition. A domain Q is also defined: 

n=\Jn e , (ii) 

e=l 

where fi does not include element boundaries. It is also useful to define the 
'interior' boundary T, 

rib 

f = U r - ( 12 ) 

i=l 

where I\ is the ith interior element boundary and is the number of internal 
inter-element boundaries. 

Consider now the function spaces S h , V h and W , 

S h = {u'i e H 1 (fi) I ttftn. G P kl (n e ) V e, Ul = ft on T g } , (13) 
V' 1 = {«# e ff 1 (O) | ti^ln. e P fcl (« e ) V e, i«i = on T g } , (14) 
W' 1 = {q h e L 2 (n) | q h \n e G P fe2 (fi e ) V e} , (15) 

where Pfc represents the space of polynomial finite element shape functions of 
order k. The spaces S h and V h represent usual, C° continuous finite element 
shape functions. Note that the space W h contains discontinuous functions. 

3.1 Standard Galerkin weak form 

The standard, continuous Galerkin problem for the equilibrium equation Q is 
of the form: given e E W h , find u h e (S h ) n such that 

/ V s w h : (1-D (e h )) C:V s u h dQ - [ w h ■ h dT = V w h G (V h ) n , (16) 
Jn Jr h 

where it was already assumed that u h is C° continuous (see equation QliSfO. A 
second Galerkin problem is constructed to solve for e (equation @). It consists 
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of: given u h E (S h ) n , find e h € W h such that 



f q h e h dQ- q h e^ d£l+ [ Vq h • c 2 Ve* q dn 
Jn Jn in 

- J gVVe^ ■ n dT = V q h eW h . (17) 

Recall that discontinuities in q h and e h across T are permitted. At this stage, the 
only boundary conditions implied by the formulation are on the displacement 
field (by construction) and the traction. Discussion regarding the enforcement 
of 'non-standard' boundary conditions is delayed until the following section. 

Two problems exist in the preceding Galerkin formulation. The first is that 
the weight function q h can be discontinuous, hence Vq h is not necessarily squarc- 
integrable on f2. This problem can be circumvented easily by requiring C° con- 
tinuity of the functions in W h . The second problem, which is less easily solved, 
is that £p q is computed from V s it'\ Therefore, calculating Ve^ everywhere in 
fi requires that the displacement field u h be C 1 continuous if singularities are 
to be avoided on T. However, since u h e H 1 (f2) (see equation ijl^ l. it is not 
necessarily C 1 continuous. 



3.2 Discontinuous Galerkin form 

The approach advocated here avoids the need for C 1 continuity of the dis- 
placement field by imposing the required degree of continuity in a weak sense. 
Before proceeding with the formulation, jump and averaging operations are de- 
fined. The jump in a field a across a surface (which is associated with a body) 
is given by: 

[a] = ai • Tii + «2 • ri2, (18) 

where the subscripts denote the side of the surface and n is the outward unit 
normal vector (ni = — n-2 in the geometrically linear case). The average of a 
field a across a surface is given by: 

ai + a 2 , 
(a) = ^ ' ( 19 ) 

Consider now the problem <)Wells et all 12004^ : given u 
W h such that 

q h e h dQ - [ q h e h C(i dSl+ I X7q h ■ c 2 Ve£ q dil - [ q h c 2 X7e^ q ■ n dT 
u Jn Jn Jr 

^M-c 2 (v £e \) dr-^(v^>-c 2 [ £ y dr 

^- lq h ] ■ [<] dT = V q h G W\ (20) 

where a is a penalty-like parameter, and h is the element dimension. No gra- 
dients of £g q or q h appear in terms integrated over Q (which includes interior 
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boundaries) in equation (|2UII . hence the continuity requirements on the spaces 
S h and W h are sufficient. 

Equation l|2U[l resembles the 'interior pen alty' method, which belongs to the 
discontinuous Galerkin family of methods l) Arnold et all l2002() . Terms have 
been added to the weak form that for a conventional elasticity problem would 
lead to a symmetric formulation. Symmetry is however not of relevance here as 
the functions q h and £^ q will generally come from different function spaces. This 
formulation is general for the case in which the space W h contains discontinuous 
functions. However, note if all functions in the space W h are C° continuous, the 
formulation is still valid, with terms relating to the jump in £g q remaining. The 
formulation would then resemble a continuous / discontinuous Galerkin method 
l|Eneel et alll2002|) . 

The solution of the gradient enhanced damage problem requires the simul- 
taneous solution of equations (|16|l and (|20|l , which are coupled. Considering the 
higher-order boundary condition (V (Ve oq ) n) ■ n = on T D Tg, the problem 
involves: find u h S (5^)™ and s h <E W h such that 

f V s w h : (1-D (e h )) C:V s tt' 1 dQ - [ w h ■ h dT - [ w h ■ f dtt 
Jn Jr h Jn 

+ a 2 h 2 (V • (V (V™") cq )) n ■ Ec 2 (V • (Ve^)) n dT 

= Vw h € (V h Y\ (21) 



/ q h e h dQ - [ q h e'* cl dQ + [ X7q h ■ c 2 Ve^ q dfl - [ q h c 2 \7e^ q ■ n dT 
Jn Jn Jn Jr 

- jT lq h j ■ c 2 (V £e \) dT - | f (V«*> • c 2 [ e JJ dT 

+ J f ^h h i-b h J dr = Vq h eW\ (22) 

where a,i is a penalty parameter which attempts to impose the boundary con- 
dition (V ( Ve e q) n) ■ n = 0, and E is Young's modulus. The enforcement of this 
boundary condition is not in the spirit of discontinuous Galerkin methods, as 
the resulting Galerkin problem is not consistent, in the same sense that penalty 
methods are not consistent (the restoration of consistency would require an ap- 
proach analogous to Nitsche's method). However, this is of no consequence for 
the considered problems due to a fortuitous choice of interpolation, as will be 
shown later. Preserving consistency while imposing the non-standard boundary 
condition weakly is complex due to the nonlinear dependency on e. This is 
a topic requiring further investigation. Boundary conditions are not explicitly 
applied on Td — T as the necessary conditions are implied implicitly in the formu- 
lation, as will be shown in examining consistency of the proposed formulation. 
The weak enforcement of the non-standard boundary condition proposed here 
has not been tested numerically. 



7 



The nonlinear equations in l(21|) and l(22|) are coupled through the dependence 
of D on e h and the dependency of e h on u h . Linearisation of these equations is 
discussed in I Wells et ail (|2004(h 

3.3 Consistency of the discontinuous formulation 

Having added non-standard terms to the weak form, it is important to prove 
consistency of the method. Applying integration by parts to the integral over 
fl in equation ((20(1 yields 

/ Vq h ■ c 2 Ve^ q dQ = - [ q h c 2 V 2 e^ dfl+ [ q h c 2 Ve>> q ■ n dT 
Jn Jn Jr 

+ j_ (q h ) c 2 [V4J dT + j_ lq h j ■ c 2 <V<) dr. (23) 

Inserting this expression into equation 12U|) , and employing standard variational 
arguments, the following Euler-Lagrange equations can be identified: 

e - e oq - c 2 V 2 e cq = in Q, (24) 

c 2 [ £eq ] • n = on F, (25) 

c 2 lVe cq ] =0 on F. (26) 

Equation (|24|l is the original problem over element interiors (see equation (JHJ). 
Equations 1)25(1 and 1(26(1 impose continuity of the corresponding fields across 
element boundaries. The Galerkin form (equation 1(20(1 ') can therefore be seen 
as the weak imposition of these Euler-Lagrange equations. 

4 Analysis and examples 

The formulation outlined in Section [3] has been analysed, and tested for various 
problems. Elements are differentiated on the basis of the interpolation order 
for u h (which is always C° continuous), the interpolation order for e , and the 
continuity of ~£ h . For example, an element with cubic shape functions for u h 
and C° continuous, quadratic shape functions for £ h is denoted P 3 /P 2 (C°). An 
element with quadratic shape functions for u h and discontinuous, linear shape 
functions for £ h is denoted P 2 /P 1 (C _1 ). 
For one-dimensional examples, e cq = s. 

4.1 Convergence for the elastic case in one dimension 

For an elastic problem, there exists a one-way coupling between the two weak 
equations ((21(1 and ((22(1 . In this section, the approximation of e, for given a 
solution to equation ((21|) . is examined. For an elastic bar, the fundamental 
problem is a second-order differential equation, and the second weak equation 
provides a projection of the solution to the equilibrium equation (and its relevant 
derivatives) onto the basis for e. 
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Figure 1: Element patch for error analysis of P 2 /P 1 (C°) formulation. 




Figure 2: Element patch for error analysis of P 2 /P 1 (C 1 ) formulation. 



4.1.1 Error analysis for the mixed strain field 

Here, an analysis for the L 2 -error in e h , \\e h — e — c 2 e :XX \\, where e is the exact 
strain field, is performed. In order to compare with numerical solutions, we 
have solved an elastic problem with a forcing function chosen such that s has 
a localised character to it. The analysis rests upon the crucial result that the 
one-dimensional finite element displacement field, u h , is nodally-exact when th e 
shape functions are derived from Lagrange polynomials l|Strang and FixL Il973j) . 

Consider the interpolation combinations P k+1 /P k (C°) and P k+1 / P k (C~ 1 ). 
Let the nodal values of e h corresponding to element e be denoted by <$fc( e -i)+i, 
. . . , Ske+i for the C° case, and by £>(fc+i)( e -i)+i, • ■ ■ , <5(fc+i)e for the C -1 case. 
The values of u h at the nodes corresponding to this element are <i(fc+i)( e -i)+i, 
. . . , <i(fc + i) e+ i. The nodal displacements coincide with the exact solution. The 
discrete matrix form of equation l)22|) reveals the stencil relating the 8- and d- 
valucs. In the C° case, {4( e -i)+i> • • • > <Wi} depend on {rf( fc+1) ( e _ 3 ) +1 , 
d(fc+i)( e +2)+i}> and in tne C" 1 case, {$ (fc+1 )( e _ 1 ) +1) 6 {k+1)e } depend on 
{d(fc+i)( e -3)+i, • • • , d(fc+i)(e+2)+i}; i- e -, on 5fc + 6 d- values. 

Consider a nodally-exact, 5(fc+l) th -order polynomial interpolant of the exact 
solution. This field is denoted by u p . Figure shows a patch of 5 elements, 
nodes for the d- and 8- degrees of freedom, u cx and u p for the P 2 /P 1 (C°) 
interpolation combination. Figure Eldepicts the situation for the P 2 /P 1 (C _1 ) 
interpolation combination. 

Let s p denote the strain field arising from vP . Using the triangle inequality 
the error can be bounded as follows: 

\\e h - e - c 2 e, xx \\n < \\e h - e p - c 2 e%\\n + ||e p + c 2 e p xx - e - c 2 e, xx \\n- (27) 



Since u p is of or der 5(k + 1), and i s nod ally-exact over fi e , finite element inter- 
polation theory iOde n and CarevUl984ft leads to the result 

He" - c 2 e p xx - e - c 2 e, xx \\ n < Ch 5 ^~ 2 , (28) 

where C is a constant, independent of h. 

The first term on the right-hand side in (|27|) is estimated using the stencil for 
each interpolation combination. The 5(fc + l) th -order polynomial interpolant, 
u p (x), can be written as 

5(fe+i) , 1 \™ 

u p {x) = Y, K[x-{e+-)h\ , (29) 

n=0 ^ ' 

where (e + l/2)h is the midpoint of O e with the first node of the mesh at x = 0. 
The nodal exactness of u h allows {dk( e -3)+i> ■ ■ ■ A(e+2)+i} to be written in 
terms of A n . On combining with the stencil, {Sk( e -i)+i, i5fc e +i} (for the 
C° case) and {<5(/c+i)( e -i)+i, ■ ■ ■ > £(fc+i)e} (for the C _1 case) can therefore be 
expressed in terms of A n . Using the interpolation functions for e , it is then a 
trivial matter to exactly evaluate \\e h — e p — c 2 e p xx \\n e - 

For convenience, \\e h — e p — c 2 e p xx \\Q was evaluated for each interpolation 
combination. The results are: 

• P 3 /P 2 (C°) 

\\e h ~e p - c 2 e%\\l e = (JLa, 2 + 2A 3 A 5 c 2 + ^A 5 2 /) h 5 + 0(h?); 

(30) 

• P 2 /P 1 (C°) 

\\e h - e p - c 2 e p xx \\l = (^A 3 2 + 22A 3 A 5 c 2 + 120A 5 2 c 4 ) h 5 + 0(h 7 ); 
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(31) 



• p 3 /p 2 (c- 1 ) 



- h !1 - < 2 -\ , lln = - + 1« 2 j M 2 ' -i> ' + (nirv. c!2) 



P 2 /P 1 (C~ 1 ) 



\e h -e p - c 2 e%\\l e = ^-A 2 2 - -A 2 A 4 c 2 + 3A 4 VJ h 3 

+ [\a 2 A a c 2 - 6^ 4 2 c 4 ] ah 3 + 3A4 2 c 4 a 2 h 3 ; and (33) 
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\\e h -e p - c 2 E%f ne = 36 (l - 2a + a 2 ) A 3 2 c 4 h + 0(h 3 ). (34) 

For each case listed above, the higher-order terms in h have a finite maximum 
power, 0(h l ), and h < m(Q), where m(f2) is a measure of the length of the total 
domain. Therefore it follows that there exist constants C\, C2, C3, C31, £4,65, 
C51, such that for the P 3 /P 2 (C°) element 

\\s h -eV-c 2 el x \\ 2 ne <G\h 5 , (35) 
for the P 2 /P 1 (C°) clement 

\\e h -eV-c 2 e%\\l e <C 2 h\ (36) 
for the P 3 /P 2 (C- 1 ) element 

||^-^-cV xx ||* <(^, if " = 5 ' 5 (37) 
,xxMS2 e ^(73^3 otherwise, v ; 

for the P 2 /P 1 (C- 1 ) clement 

\\e h - e» - c 2 e%\\ 2 Qe < C 4 h 3 V a, (38) 
and for the P^P"^ 1 ) element, 

\\^-e*-c 2 e* xx \^<\fl l \ a = l (39) 
1651/1 otherwise, 

where the constants Ci, . . . , C51 are independent of h (if u is sufficiently regular) 
and clement number e (since the constants can be chosen to be the maximum 
over all elements). For elements of a uniform size, in each mesh there are 
n c \ = m(Sl) jh elements. Therefore, 



11^ - ^ P - <?e%\\l < ^r 1 ™*\\e h - £ p - c 2 e%\\ 2 nc , (40) 



m(0) 
T 

leading to, for the P 3 /P 2 {C°) element 

He* - £ p - c 2 e%\\l < m(Q)Ci/i 4 , (41) 
for the P 2 /P 1 (C°) element 

||^- £ P-cV xx ||^<m(0)C 2 /i 4 , (42) 
for the P 3 /P 2 (C 1 ) element 

^- £p - c hji^<(1S~ 3/i 1 if t r 5 ' 5 ^ 

\m(\l)C3ih otherwise, 
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Element type 



P L /P»(C- L ) {a = 1) 
P 1 /P°{C- 1 ) (a^l) 

p2 / pl (C 0) 

P 2 /P 1 (C- 1 ) 
P 3 /P 2 (C°) 

p3 /p2(c -i ) ( a = 5 ,5) 
P 3 /P 2 (C- 1 ) (a =£5.5) 



convergence rate 



TJ7F) 
0(/i ) 

O /i 2 



Table 1: Analytical convergence rates in terms of ||e — e^Hn for various elements. 



for the P 2 /P 1 (C- 1 ) element 

||^_ £ P_ c 2 £ pj|2 e < m(0) g 4?l 2 

and for the P 1 /- P °(C _1 ) element 



1 "i(fi)(75i/i otherwise. 



(44) 



(45) 



Since k > in (|28[) . each of equations l|41|l- 145|l can be combined with l)28[l. 
and constants Ci, C2, C3, C31, C4, C5, C51 can be found, independent of h, 
such that for the P 3 /P 2 (C°) element 



||^_ £ P_ C 2 £ PJ|2 ^^d/i 4 , 



for the P 2 /P 1 (C°) element 



for the P 3 /P 2 (C- V ) element 



||^- £ P-cV x J|2 < 



I m 
I m 



(ft)C 3 /i 4 if a = 5.5 
(f2)C3ift. 2 otherwise, 



for the P 2 /P 1 (C- 1 ) element 

\\e h -s p -c 2 e% x f n <m(il)C 4 h 2 V a, 
and for the pi/P^C" 1 ) element, 

\m(tt)C 5 h 2 if a = 1 



' \m(Cl)C5ih° otherwise. 



(46) 
(47) 

(48) 

(49) 

(50) 



The results of the analysis are summarised in Tabled 
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Figure 3: Quadratically tapering bar. 



4.1.2 Observed convergence rates 

The convergence rate estimates from the previous section are now compared with 
observed rates. To do this, the convergence of e is examined for the quadratically 
tapering bar shown in Figure |3J Considering that the origin is located as the 
centre of the bar, the cross section A is given by 

A(x)=A 1+1 2 A 1 ^, (51) 

where 7 controls the ratio between At and A2 ('f 2 = (A2 — At) /At). This 
problem can be cquivalcntly formulated as a rod with unit cross-section and 
loading / which varies along the rod, which allows the convergence analysis 
from the previous section to be carried over for this case. The exact solution 
for e along the bar is: 

- PL2 
£ ~ EAt {L 2 + j 2 x 2 ) ' ( ' 

and the exact solution of e xx is: 



2PL 2 7 2 (37-z 



2„2 _ T? .\ 

(53) 



EAt{L 2 +j 2 x 2 y 

The term e h involves both e and s >xx . The s component involves the standard 
L 2 projection of e h (coming from the solution of the equilibrium equation) onto 
the basis for e h . This projection does not involve any of the element interface 
terms that have arisen in the discontinuous Galerkin formulation. Of special 
interest is the approximation of the e tXX term. To examine numerically the 
convergence of this term, the weak form corresponding to: 

e-e,xx = (54) 

is solved, as it is the e jXX term that dominates the convergence rate. However, 
in practice, the convergence rate of e may appear to dominate due to a large 
difference in the constants in the error inequality. Taking equation (|22|l and 



13 




/ q h e h dQ + [ q%e% dSl - [ q h e%n dT - [ \q h \ (e%) dT 
Jn Jn Jr Jr 

- jf (q%) ¥ h l dr + J f i M PI dT = v i h e ( 5fi ) 

For the parameters Ai = 1 mm 2 , A-i =0.1 mm 2 , L = 50 mm, E = 1 MPa 
and P = 1 N, the convergence behaviour is examined for a range of elements. To 
gain insight, the form of the exact solution is shown in Figure^J The L 2 -norm of 
the computed error for a range elements is shown in Figure[S]for the case a = 1. 
If the error is computed by integrating over the entire bar (—50 < x < 50), 
as in Figure |SJ the results are polluted by errors at the boundaries of the bar. 
In Figure the convergence behaviour represents primarily how rapidly the 
computed solution approaches the exact solution at the boundaries. However, 
when simulating a tapered damaging bar, the error at the boundaries of the 
bar is of little consequence as damage develops at the centre, and e at the 
boundaries plays no role (presuming that the deviation from the exact solution 
is not sufficient to induce spurious damage development). Excluding the error 
at the end of the bar by integrating over —40 < x < 40, the convergence 
behaviour is more predictable, as can be seen in Figure [S] From Figure it 
can be concluded that the convergence rate for all elements is consistent with 
the predicted rates (as summarised in Tabic QJ. 

The effect of a on convergence for the P 2 /P 1 (C _1 ) element is shown in 
Figure [7| As expected, the convergence rate is unaffected by a. For the 
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Figure 5: Convergence for the problem e 
-50 < x < 50. 



with a — 1 on the domain 




Figure 6: Convergence for the problem e = s yXX with a = 1 on the domain 
-40 < x < 40. 
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Figure 7: Convergence of the P 2 / P 1 (C x ) element for different values of a on 
the domain —40 < x < 40 for the problem e = e„. 



P 3 / p2 (C^ 1 ) element, the convergence behaviour for different a is shown in 
Figure |H1 From Figure |H1 it is clear the convergence rate increases by one order 
for a — 6. This is close to the predicted value of a — 5.5. 



4.2 Convergence for the inelastic case in one dimension 

The formulation is now examined for the inelastic case. Again, the quadratically 
tapering bar is considered, which leads to damage development at the centre of 
the bar. For a particular relationship between D and n, it is possible to solve 
the problem analytically, which provides the basis for numerical convergence 
tests. 



4.2.1 Analytical solution 

Consider again the bar shown in Figure where the shaded region indicates 
the damaged zone and x = ±w is the location of the damage-elastic boundary. 
For convenience, the force P is expressed as 

P = (1 - I3 2 ) EA lKo , (57) 

where f3 governs the magnitude of the applied load and no is the value of K at 
which damage is first induced. In the undamaged part of the bar, the strain 
response is given by 

£ = 4 ^ 

and within the damaged zone by, 

e =(rdbz- (59) 
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Figure 8: Convergence of the P 3 / P 2 (C x ) element for different values of a on 
the domain —40 < x < 40 for the problem e = s xx . 



With the intention of finding an analytical solution, a simple damage law is 
assumed, 

{0 if K < Kn 

, «o ., ^ (60) 
1 it K > Kq. 
K 

The relationship in equation (|60[) is the analogy of perfect plasticity for the case 
c = 0, in the sense that it yields a plateau in the load-displacement response 
once the elastic limit has been exceeded. Assuming that no unloading takes 
place in the damaged zone, k = e, that is 

k = e + c 2 e. xx — w < x < w. (61) 

Inserting equations (|(jlfl and H(jU|) into equation (|59|) . the following ordinary 
differential equation is obtained: 

P \ P 
1 — )e-c 2 — — e, xx =0 -w<x<w, (62) 

which governs the response in the damaged zone. The general solution to the 
above equation for x > is given by: 

M( w ( ?L I 



Uc 7N /i+F'4' Lcy/T+p> J Uc 7N /i+F' 4' Lc-y/T+fP 
£ = Cl ^ + ° 2 ^ ' (63) 

where M and W are Whittaker's functions l Abramowitz and StegurJ . ll96. t T() . and 
C\ and Ci are integration constants. The integration constants can be obtained 
by considering boundary conditions at x = ±w. The first considered condition 
is symmetry about x — 0. Expanding equation l|63|) in a Taylor series about 
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Figure 9: Error in e h for the damaging bar. 



x = 0, and requiring that odd terms vanish, the following relation is obtained 

2tt 



C\ — C2 



p f 1 _ -P 2 L 2 + C1 ^l+!3 2 \ ' 
4 4cVl+7/3 2 



(64) 



where T is the Gamma function. The second considered condition is: 

[e J =0 atx = w, (65) 

which implies 



£ ' x ~ \EA 



at x = w. 



(66) 



Finally, the location of the elastic-damage boundary can be determined by 
requiring that 

k = kq at x = w. (67) 
For the sake of brevity, the expressions for C\, Ci and w have been omitted. 

4.2.2 Damage convergence results 

The convergence behaviour of various elements for the outlined test problem 
is now examined. For the tests, the following parameters are adopted: 7 = 3, 
(3 = 5/29, L = 100 mm, E — 200 x 10 3 MPa, c = 1 mm and Ax = 1 mm 2 
and kq = 1 x 10~ 3 . In calculating the error, the Gamma function has been 
computed numerically. 

The results of the error analysis, for various elements, are shown in Figure El 
For the P 1 /P° (C*- 1 ), P 2 /P 1 (C^ 1 ) and P 3 /P 2 (C^ 1 ) elements, a = 1, a = 4 
and a = 6, respectively. These choices draw upon the observed convergence 
results for the elastic case. For all the elements, the solution converges. Once a 
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level of refinement has been reached, the convergence rate is linear. Moreover, 
for a given interpolation order, continuous and discontinuous interpolations yield 
similar results. 



4.2.3 Non-trivial damage response 



The quadratically tapering bar is now examined for a non-trivial softening rela- 
tionship. For the tapered bar the relevant parameters are: 7 = 3, A\ — 1 mm 2 , 



Young's modulus E = 20 x 10 3 MPa, n = 1 x 10" 



= 0.0125, and c = 1 



The functional form of the damage variable is specified to be 



D=ll- 



Kq(k c - K) 

k(k c - k ) 



if K < Kq 

if K < K < K c 

if K > K c , 



(68) 



which leads to a linear softening relationship for c = 0. The nu merical per- 
forma nce of the formulation has previously been demonstrated in IWells et alJ 
l)2004[) for continuous, piecewise linear u h and constant e h , discontinuous across 
element boundaries (the P 1 /P°(C _1 ) element). The formulation was shown 
numerically to converge to a benchmark solution. Here, the formulation is ex- 
tended for a range of different element types. 

The motivation behind the considered strain gradient-dependent model is 
regularisation in the presence of strain softening. Without strain gradient effects 
(c = 0), computed results are pathologically mesh-dependent; the result of which 
is manifest in the load-displacement responses. Therefore, each of the elements 
which to this point have been examined are tested, and the load-displacement 
responses reported for meshes with 20, 40, 80, 160 and 320 elements. For the 
P 3 /P 2 (C _1 ) element, the numerical tests are performed with meshes of 20, 40, 
80 and 160 elements. To provide a reference solution, the response computed 
using 160 P 3 /P 2 (C°) elements is included in all figures. 

The load-displacement responses for the two elements using a continuous in- 
terpolation of e are shown in Figure [Till As the mesh is refined, the computed 
response for both element types converges towards to the reference solution. 
The load-displacement responses for three elements using a discontinuous inter- 
polation of e are shown in Figure EU Again, for the P 1 / P° (C" 1 ) element, 
a = 1, for the P 2 /P 1 (C^ 1 ) clement, a = 4, and for the P 3 /P 2 (C^ 1 ) element, 
a = 6. It is clear, for all elements, that the load-displacement response con- 
verges to the reference solution with mesh refinement. To further examine the 
computed results, the damage profiles for the two continuous elements and the 
three discontinuous elements are compared in Figures 1121 and 1131 respectively. 
For all cases, the 160 element mesh is considered. Clearly, the damage profiles 
are nearly identical for all element types. 
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Figure 12: Damage profiles for the e^-continuous elements, computed using 160 
elements. 
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Figure 13: Damage profiles for the e^-discontinuous elements, computed using 
160 elements. 
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Figure 14: Three-point bending specimen. 



4.3 Three-point bending test 

To conclude the numerical validations, a three-point bending test is performed 
using a P 2 /P 1 (C°) element. The element is triangular, with degrees of freedom 
for u located at the vertexes and at the mid-sides, and degrees of freedom 
for e h located only at the vertexes of the element. The choice of a quadratic 
interpolation of u h means that the penalty term for imposing the non-standard 
boundary condition in equation (|21|l vanishes, and the non-standard boundary 
condition is satisfied by construction (see equation 11211) ). The equivalent strain 
is taken as the trace of the strain tensor, 

e oq = trace (e) . (69) 

This choice does not reflect a strong physical motivation, rather it is chosen for 
illustrative purposes as it allows for relatively simple linearisation of the method 
(which can become extremely complex when c =/= 0). 

The three-point bending test is performed for two different meshes with c ^ 
and c = 0. The adopted geometry for the three-point bending test is shown in 
Figure El an d the adopted material parameters are: Young's modulus E — 
20 x 10 4 MPa, Poisson's ratio v = Q, k = 1x 1CT 4 , and k c = 1.25 x 1CT 2 . For 
gradient-dependent simulations, c = 8 x 10 -2 mm. Computations are stopped 
when damage reaches unity at any point in the mesh, and damage development 
at the supports is prevented. The computed damage contours for the four cases 
are shown in FigureEl From the damage contours, it is clear that the computed 
results are similar for the two meshes with c ^ 0. In the absence of regularising 
effects (c = 0), the result is clearly affected by the discretisation. The load- 
displacement responses for the various cases are shown in Figure El Recall 
that a computation is halted when damage reaches unity at a material point. 
For the gradient-dependent case, the responses for the two meshes are similar. 
For the case c = 0, the responses are also similar, which is somewhat in contrast 
to what is normally expected for a strain softening problem. The responses are 
similar in this case due to the spurious development of two cracks for the finer 
mesh (in contrast to the single main crack for the coarse mesh). This is evident 
from the damage contours in Figure El 
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(a) (b) 

Figure 15: Damage contours for two meshes (a) with gradient effects (c 7^ 0) 
and (b) without gradient effects (c = 0). 
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5 Conclusions 



A discontinuous Galerkin formulation for a strain gradient-dependent damage 
model has been investigated for a range of different finite elements. The model 
allows the numerical solution of a continuum problem which would classically re- 
quire C 1 interpolations with a simple C° or even discontinuous basis. Examples 
demonstrate robust performance for a range of polynomial orders and degrees 
of continuity of the interpolation functions, and are supported by rigorous error 
analysis. Specifically lower-order interpolations perform well and are relatively 
simple to construct. The convergence properties of the proposed method have 
been examined for the elastic case, for which the observed rates are consis- 
tent with the theoretically predicted rates. The formulation has been observed 
numerically to converge also for damage problems. Finally the formulation 
was applied successfully to a two-dimensional problem. While the approach is 
promising, several issues remain. Difficulties which must be resolved for other 
gradient models include the effective imposition of boundary conditions on the 
fixed boundary of a body, and at moving boundaries internal to a body. The 
development of thermodynamically consistent models would assist in this sense, 
as the higher-order kinematic gradients have a natural partner in the energetic 
sense. 
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